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^ ■ ABSTRACT 
> ■ 

O ■ Integrations in fixed N-body realisations of smooth density distributions corresponding to 

a chaotic galactic potential can be used to derive reliable estimates of the largest (finite time) 
! Lyapunov exponent xs associated with an orbit in the smooth potential generated from the 



(N 



same initial condition, even though the A^-body orbit is typically characterised by an A-body 
J> ' exponent XN ^ XS- This can be accomplished either by comparing initially nearby orbits in 

■ a single A-body system or by tracking orbits with the same initial condition evolved in two 

, different A-body realisations of the same smooth density. 

I> 

o 

^ 1. INTRODUCTION AND MOTIVATION 

^p^i At the present time, there are two general approaches which can be used to model the 

O ' structure and evolution of systems like elliptical galaxies. On the one hand, one can perform 

' detailed A-body simulations, solving the coupled equations of motion either exactly or in some 

Ci • approximation. On the other, one can try to construct equilibrium models as solutions to the 

collisionless Boltzmann equation, either analytically (e.g., using the approach developed by 
Hunter & Qian [1993]) or numerically {e.g., by implementing some version of Schwarzschild's 



X 



; [1979] method). 

^ Both approaches would appear to be extremely useful, but each has its limitations. Despite 

significant advances in hardware, direct A-body simulations still cannot be performed for par- 
ticle number as large as A ~ 10^^, so that one cannot consider point masses in the integrations 
as corresponding to individual stars in real galaxies. Integrations in smooth potentials have the 
advantage that one can consider characteristics which, presumably, correspond to the orbits of 
individual stars. However, the very assumption that the system can be described by a smooth 
potential constitutes an idealisation which, albeit generally accepted, has never been proven 
rigorously. 

Obvious questions to be answered thus include the following: To what extent is it true that, 
for sufficiently large A, solutions to the full gravitational A-body problem can be mimicked by 
motions in a smooth potential? And, to the extent that the smooth potential approximation is 
not completely sufficient, to what extent can discreteness effects really be modeled by friction 
and noise in the context of a Fokker-Planck description (c/. Rosenbluth, MacDonald, & Judd 
1957)? In particular, to what extent are 'real' A-body orbits well mimicked by solutions to a 
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(time-dependent) Langcvin equation {cf. Chandrasckhar 1943) which incorporates dynamical 
friction and Gaussian wliitc noise? Fokkcr-Planck descriptions were formulated originally to 
extract statistical information about long time behaviour, assuming implicitly that the bulk 
potential is integrable or nearly integrable. However, recent years have seen a growing recogni- 
tion that galactic potentials may admit a fair amount of chaos; and analyses of the short time 
behaviour of individual orbits in chaotic potentials have provided compelling evidence that 
friction and noise can dramatically accelerate phase space transport (Licbcrman & Lichtcn- 
berg 1972, Lichtenberg & Wood 1989, Kandrup, Pogorelov, &; Sideris 2000, Siopis &: Kandrup 
2000). Does this really mean that discreteness effects can be important already on time scales 
much shorter than the relaxation time Ir"! 

Closely related to these issues is the nature of the continuum limit. In what sense is it 
true that, as ^ oo, orbits in the A^-body potential converge towards characteristics in some 
smooth potential? Superficially, at least, it might seem that such a convergence is impossible. 
Dating back to Miller (1964), it has been recognised that the AT-body problem is chaotic in 
the sense that individual orbits exhibit exponential sensitivity towards small changes in initial 
conditions; and it seems generally accepted today that, when expressed in units of inverse 
dynamical times i^^, the largest A^-body Lyapunov exponent xn does not converge towards 
zero as N ^ oo {cf. Kandrup & Smith 1991, Goodman, Heggie &; Hut 1994), even if the N- 
body system samples an integrable density distribution. Indeed, recent work, both numerical 
(Hemsendorf & Merritt 2002) and analytic (Pogorelov 2001), suggests that, even for a density 
distribution corresponding to an integrable potential, the largest Lyapunov exponent may 
actually increase with increasing A^. As probed in the usual way, the A"-body problem may 
become more chaotic as N increases! 

A complete resolution to this apparent conundrum will require long time integrations of 
systems with very large A^, which is impractical using current hardware. However, considerable 
insight into the continuum limit can be, and has been, obtained by studying the properties 
of orbits and orbit ensembles evolved in frozen-N systems, i.e., fixed (in time and space) AT- 
body realisations of specified smooth density distributions (Kandrup &; Sideris 2001, Sideris & 
Kandrup 2002). In particular, that work led to several significant conclusions: (1) The largest 
A^-body Lyapunov exponent xn does not decrease with increasing A^, even for an integrable 
density distribution. However, there is still a clear, quantifiable sense in which, as A^ increases, 
the AT-body orbits become progressively more similar to smooth potential characteristics. (2) 
As N increases the Fourier spectra associated with AT-body orbits more closely resemble the 
spectra associated with characteristics in the smooth potential, be these either regular or 
chaotic. (3) Alternatively, viewed macroscopically A^-body orbits and smooth characteristics 
with the same initial condition typically diverge as a power law in time on a time scale tciN) 
that increases with increasing N. For the case of regular characteristics, to oc N^/'^to', for 
chaotic characteristics, to (InAQt/j. It follows that, for sufficiently large N, A^-body orbits 
and smooth potential characteristics remain close for comparatively long times. 

This would seem a result of some significance, but it still begs an important issue. If A''- 
body orbits converge towards smooth potential characteristics, it should be possible, at least 
for sufficiently large to extract information about any chaos that may be associated with 
the bulk potential. In particular, it must be possible to extract estimates of finite time {cf. 
Grassberger, Badii, & Politi 1988) Lyapunov exponents xs for motion in the smooth potential, 
even though the A^-body orbits themselves are characterised by exponents xn which, typically, 
are much larger than xs- (Chaotic orbits in generic smooth potentials typically have a largest 
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Lyapunov exponent xs ~ tjj^ . For both interacting [Hemsendorf & Mcrritt 2002] and frozen- iV 
[Kandrup &: Sidcris 2000] Plummer systems, XN ~ 20t^^ for iV ~ 10^.) 

The aim of this note is to demonstrate that this can in fact be done. In particular, it 
will be shown that there are at least two different ways in which estimates of xs can be 
extracted from frozen- A?^ systems, one involving a comparison of orbits generated from nearby 
initial conditions and the other involving integrations of the same initial condition in two 
different A^-body realisations of the same smooth density distribution. Section 2 describes the 
algorithms and then exhibits representative results for two simple model potentials. Section 3 
interprets the success of these algorithms by postulating the existence of two 'types' of chaos, 
microscopic chaos, or microchaos, which is generic to the A?^-body problem, and macroscopic 
chaos, or macrochaos, possibly associated with the bulk potential which, if present, will also 
be manifested in the A'^-body problem. 

2. NUMERICAL ESTIMATION OF SMOOTH POTENTIAL LYAPUNOV EXPONENTS 

2.1. Models Considered 

The algorithms described in this Section will be applied to representative orbits evolved in 
two different density distributions: 

1. A spherically symmetric Plummer sphere, for which 

This corresponds via Poisson's equation to a spherically symmetric, and hence integrable, 
potential 

GM 

Units were so chosen that G = M = b = 1. 

2. A constant density triaxial ellipsoid, for which 

, , 3M fm^ if m2 < 1. 
^^^^) = 4;^n0 ^im^>l, 

with 

(4) 




perturbed by a spherically symmetric density spike (black hole) of mass Mbh , this correspond- 
ing to a potential 

^E{r) = ^o + l {u^y + ulv^ + ulz^) - (5) 

with e = 10-3. Attention focused on the case M = 1.0 and Mbh = lO-^-^M k, 0.0316228, 
and units were again chosen so that G = 1. The axis ratios were taken as a = 1.95, h = 1.50, 
and c = 1.05, which yield (c/. Bertin 2000) $o ~ - 1.00608, Ua ~ 0.4663, ojb ~ 0.5508, and 
ujc ~ 0.6753. For energies sufficiently small that orbits are restricted to m < 1, the phase 
space is almost completely chaotic (Kandrup &: Sideris 2002). This implies that one need 
not worry about transitions between regular and chaotic behaviour which can be induced by 
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discreteness effects in more complex potentials like the triaxial generalisations of the Dehnen 
(1993) potential, which admit a complex coexistence of both regular and chaotic orbits. 

Frozen- A'^ orbits in these systems were integrated with a particle-particle numerical scheme 
using a variable timestep integrator with accuracy parameter 10~^ which conserved energy to 
at least one part in 10^. The 1/r kernels for the individual masses were regulated through the 
introduction of a softening parameter e = 10""^. 

3.2. The Numerical Algorithms 

Algorithm 1. Involving a comparison of orbits generated from nearby initial conditions in the 
same frozen-iV density distribution. 

The key realisation here is that, even though orbits generated from nearby initial conditions 
diverge initially at a rate A ~ xn^ tliis divergence quickly saturates once the growing separation 
between the orbits becomes large compared with the typical distance between neighboring point 
masses. If the smooth potential characteristic is regular, this initial exponential divergence is 
replaced immediately by a more modest power law divergence which proceeds on a time scale 
tc oc N^^'^t£). If instead the characteristic is chaotic, the initial exponential diverg ence at a rate 
A ~ xn is replaced by a slower exponential divergence at a rate A ~ X5 which, in turn, typically 
proceeds until the separation becomes 'macroscopic', i.e., comparable to the size of the entire 
accessible phase space region. At this point, this second exponential divergence saturates and 
is replaced by a power law divergence which proceeds on a time scale to oc {lnN)tD- 

These statements can be corroborated straightforwardly by tracking the separation 

,5r(t) = |ri-r2| (6) 

and plotting ln5r(t) as a function of t. By so doing, one discovers that, just as for the 
smooth potential, nearby initial conditions can exhibit significantly different values of xs- 
(Nearby initial conditions evolved in a smooth potential can have significantly different finite 
time Lyapunov exponents even if, for late times, these exponents converge towards the same 
asymptotic Xoo-) If) alternatively, one wishes to estimate a 'typical' xs for orbits in some given 
phase space region, one can select an ensemble of initial conditions from that region, evolve 
each initial condition into the future, and then extract a mean {xs) from the time-dependent 
mean separation 

k 

W = ^E^^- (7) 

i=l 

Figure 1 exhibits the results of such a computation for an ensemble of 100 chaotic initial 
conditions evolved in a frozen- realisation of the the ellipsoid plus black hole system. The 
different curves in the Figure represent frozen- A?^ backgrounds with N varying between N = 
10^-^ and N = 10^. Figure 2 exhibits analogous data generated for regular initial conditions in 
the Plummer potential. In each case, the initial conditions sampled a phase space region of size 
Ar ~ Av ~ 10^'^, with the perturbed orbits being generated from initial conditions displaced 
from the original initial conditions by a distance 6r{Q) = 10~^ in a randomly chosen direction. 

It is evident that, for both the ellipsoid and Plummer potentials, the mean separation (Sr) 
begins by diverging at a rate that is comparable to the value of the largest A?^-body Lyapunov 
exponent, XN- For the case of the ellipsoid potential, this exponential divergence is (at least 
for sufficiently large N) eventually replaced by a slower exponential divergence at a rate ~ xs 
which persists until 6r becomes 'macroscopic,' i.e., comparable to the size of the accessible 
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phase space region. At this point the exponential divergence is replaced by a slower power 
law separation which proceeds until 5r saturates. For the case of the Phimmer potential, the 
second exponential phase is absent, the initial exponential phase being replaced immediately 
by a power law growth. This is especially evident from Figure 3, which plots the same data as 
Figure 2 on a linear scale. 

As noted already, for both regular and chaotic potentials the first exponential phase 
typically stops once the separation 5r becomes large compared with the mean interparti- 
cle spacing ~ n~^^^, with n a typical number density. For the constant density ellipsoid, 
n-V3 ~ 0.427iV-V3. por N = 10^ this corresponds to n^Vs ~ 0.00919 and Inn^Vs ~ -4.69; 
for N = 10~^, n~^/^ 0.00427 and lnn~^/^ — 5.45. Alternatively, the second, slower expo- 
nential phase exhibited by chaotic potentials typically ceases once the separation has become 
'macroscopic', i.e., comparable to the size of the accessible phase space region. This macro- 
scopic scale appears to coincide with the scale on which two initially nearby orbits evolved in 
the smooth potential will cease their exponential divergence. 

As discussed more carefully in Section 3, the comparatively sharp break between the first 
and second phases indicates that one can implement a de facto distinction between microscopic 
chaos associated primarily with close encounters and macroscopic chaos associated with the 
bulk potential, even though the effects of these sources of chaos are not completely decoupled. 

It is evident that this prescription for estimating xs can only work for comparatively large 
values of N, where the typical interparticle spacing is much smaller than the total size of 
the accessible configuration space region. If this condition is not satisfied, 5r{t) will become 
'macroscopic' almost as soon as it becomes large compared with n~^/^, so that the intermediate 
second stage disappears. For the particular model exhibited in Figure 1, one requires N > 10^ 
or so in order to obtain a reasonable estimate of xs- Indeed, it is evident from the dot-dashed 
curve in Figure 1 that, after the initial exponential divergence and before the final saturation, 
the data for N = 10^'^ can be well fit by a linear growth law 5r = A{t — to). 

Algorithm 2. Involving a comparison of orbits generated from a single initial condition evolved 
in two different frozen- density distributions that sample the same smooth density. 

The smooth exponent xs also provides information about the rate of divergence associated 
with orbits generated in two different frozen- AT simulations. Specifically, if one computes 5r for 
a collection of orbits evolved in two different frozen- A^ density distributions associated with the 
same chaotic potential, one again finds an evolution manifesting the same three stages. This 
is illustrated in Figure 4, which was generated for the ellipsoid potential for the same initial 
conditions as Figure 1. Figure 5 exhibits analogous data generated for the Plummer potential. 
It is evident once again that, for the Plummer potential, the intermediate stage is absent. 

Given that xs provides information about the rate of divergence of orbits in two different 
frozen-A" backgrounds, each of which can be viewed intuitively as a 'perturbation' of the 
smooth density distribution, one might also expect that xs provides information about the 
rate at which orbits in a single frozen- AT simulation diverge from smooth characteristics with 
the same initial condition. As illustrated in Figure 6, this expectation is in fact correct. 
The macroscopic power law divergence between the frozen-A orbit and the corresponding 
smooth potential characteristic, which parallels the behaviour observed if Algorithm 1 or 2 
is implemented, has been discussed extensively elsewhere (c/. Figure 8 in Kandrup & Sideris 
2001 and Figure 2 in Sideris & Kandrup 2002 and the accompanying discussion). 
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3.3. Why these algorithms work 



Given the assumption that orbits in an AT-body system 'feel' two different sorts of chaos 
which act on different scales, it is not surprising that one can derive estimates of both the 
A/^-body Lyapunov exponent XN and the smooth potential xs from a comparison of initially 
proximate orbits in a single frozen- system. However, the fact that xs is also related to the 
rate of divergence of orbits in different frozen-iV simulations is, perhaps, less obvious. The key 
to understanding this phenomenon is the fact that discreteness effects really can be mimicked 
by dynamical friction and Gaussian white noise in the context of a Langevin description. 

Specifically (Sideris and Kandrup 2002), at both the level of individual orbits, as probed, 
e.g., by Fourier spectra, and at the level of orbit ensembles, as probed, e.g., by the efficiency 
of phase mixing for both regular and chaotic orbit ensembles, discreteness effects associated 
with an A'^-body density distribution arc extremely well reproduced by Gaussian white noise 
with a 'temperature' ~ \E\, with E the orbital energy, and a coefficient of dynamical friction 
r] oc 1/N. This dependence on A^ is of course very similar to the scaling tn ^ r]~^ oc N/{lnN) 
predicted in a conventional Fokker-Planck description; and, indeed, given the limited range in 
N that can be probed numerically (the notion of a smooth potential appears to break down 
for N <10^; simulations with A^ » 10^ become prohibitively expensive computationally!), the 
numerical simulations are completely consistent with this scaling. 

But why does this explain the observed divergence of orbits in different frozen- A" distri- 
butions? The crucial point here, as described, e.g., in Habib, Kandrup &; Mahon (1997) or 
Kandrup and Novotny (2002), is that, viewed mesoscopically, an ensemble of noisy orbits with 
fixed O and rj, each generated from the same chaotic initial condition or from a set of very 
nearby initial conditions, will typically disperse in such as a fashion that 

<5roc (e7?)V2exp(x50- (8) 
Given the assumed scaling 77 oc 1/N, it then follows that 

In Sr = const -|- ^ In 77 -|- xst = const — ^ In AT -|- xst- (9) 

Numerical simulations demonstrate that noisy orbits diverge at the same rate xs observed 
for orbits in frozen- A^ simulations; and the connection between N and r] implicit in a Fokker- 
Planck description makes a specific prediction as to the A/"-dependence of the exponential 
prefactor. If discreteness effects really can be modeled as Gaussian noise, 5r should satisfy 
eq. (9). To the extent that the 'mean' trajectory associated with the noisy ensemble coincides, 
at least approximately, with the smooth potential characteristic, the same scaling should also 
be observed when comparing noisy orbits and smooth potential characteristics. 

This scaling implies (a) that \n6r should grow linearly at a rate xs, independent of N, 
but (b) that, for fixed t, an increase in A^ or a decrease in rj by an order of magnitude should 
decrease \n6r by (1/2) In 10 ~ 1.15. That this scaling is in fact realised for both noisy orbits 
and orbits in frozen-A^ backgrounds is evident from Figure 7. Here the solid curves exhibit 
the results of noisy integrations of the same initial conditions used to generate Figure 1 and 
4, all computed with @ = 1.0 but allowing for values of 77 extending from 10~^ to 10~^. 
The dots accompanying the upper two curves represent results from frozen- A?^ integrations for 
A'^ = 10^'^ and 10^'^ (the same data plotted in Figure 6). The results from Sideris and Kandrup 
(2002) suggest a best fit correspondence log^g?? = — logj^Q A^ +p, with p « 0.5; and indeed, it 
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is apparent visually that the noisy and frozen-A/" curves for, e.g., rj = 10~^ and N = 10^-^ are 
extremely similar. 

Presuming that this scaling holds for smaller r] as well, the lowest curve in Figure 7 should 
correspond, at least approximately, to N = 10^-^. 

III. DISCUSSION 

This paper has described two algorithms which can be used to obtain estimates of the 
largest (finite time) Lyapunov exponents xs associated with a smooth density distribution 
from frozen-A/" realisations of that density. The first involves comparing two orbits in a single 
frozen- A?^ system generated from nearby initial conditions. The second involves comparing 
orbits with the same initial condition evolved in two different frozen- A^ systems, each sampling 
the same smooth density distribution. The success of the first algorithm emphasises the fact 
that detailed information about the bulk potential really is buried in an A^^-body simulation. 
The success of the latter emphasises another important point, namely that smooth potential 
Lyapunov exponents xs also provide information about the divergence of the same initial 
condition in different A^-body systems, i.e., information about the extent to which, viewed 
mesoscopically, discreteness effects limit the intrinsic reliability of orbits in a pointwise sense. 
As stressed already, the success of this alternative algorithm reflects the fact that discreteness 
effects really can be well mimicked by Gaussian white noise in the context of a Fokker-Planck 
description. 

The key point in all this is that, for the case of a chaotic bulk potential, two nearby initial 
conditions will, when evolved into the future, exhibit a three stage evolution, reflecting the 
effects of both microscopic and macroscopic chaos, i.e., microchaos and macrochaos. 

(1) For early times and small separations, the orbits will diverge exponentially at a rate com- 
parable to a typical A^-body Lyapunov exponent XN- 

(2) However, once the separation between the orbits becomes large compared with the typical 
interparticle spacing, this divergence ceases and is replaced by a slower divergence at a rate 
~ Xs, which proceeds until the separation becomes macroscopic. 

(3) At still later times the orbits exhibit a more modest power law divergence. 

The computations described here only yielded estimates of finite time Lyapunov exponents, 
not the true Lyapunov exponent as defined in a late time limit, which, given a complex phase 
space, can be much larger or smaller. This, however, is not necessarily bad. Although old in 
physical time, galaxies are young objects when expressed in terms of the dynamical time to 
- typically no more than ~ 100 — 200t/) in age - so that such asympotic limits are not well 
motivated physically. However, one might argue that, even though an asymptotic limit is not 
justified for individual orbits, the true Lyapunov exponent is important in that (c/. Kandrup 
Sz Mahon 1994) it characterises the average instability associated with the invariant measure, 
i.e., a uniform sampling of the chaotic portions of the constant energy hypersurface. The 
obvious point, then, is that to obtain an estimate of the true Lyapunov exponent it suffices 
to repeat the calculations described here for an ensemble of initial conditions sampling the 
constant energy hypersurface which, as described elsewhere (c/. Kandrup, Sideris, &; Bohn 
2001) is straightforward numerically. Alternatively, one can actually compute estimates of xs 
in the usual way (Benettin, Galgani, and Strelcyn 1976) by tracking the evolution of a small 
perturbation which is periodically renormalised, provided only that one makes sure that the 
perturbation always remains large compared with the scale on which the microscopic chaos 
saturates, i.e., very large compared with a typical interparticle spacing. 
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As a practical matter, it would appear that the algorithm can work for any A'^-body system 
in which the typical intcrparticle spacing is sufficiently small compared to the size of the system. 
If N is not sufficiently large, the first exponential phase will not saturate until the separation of 
the originally proximate orbits has become macroscopic, so that the second exponential phase 
is lost. For the models considered here, one requires N > 10^ or so. For systems manifesting 
very high density contrasts, e.g., triaxial generalisations of the cuspy Dchnen potentials, one 
may require much larger N to obtain an adequate sampling of the central region. In point 
of fact, however, the requirement of large N is more than a practical consideration: it would 
appear that, if N is too small, the very notion of a bulk potential becomes suspect. One finds, 
e.g., that, for the ellipsoid plus black hole potential, discreteness effects can be reasonably 
well modeled by Gaussian white noise for N ~ 10^, but that this model fails for substantially 
smaller N (Sideris & Kandrup 2002). 

Viewed from the standpoint of nonlinear dynamics, the A-body problem - or at least the 
frozen- A model considered here - constitutes an interesting example of a system in which chaos 
can arise for different reasons on different scales. Because gravity is strong on short scales, one 
finds generically that close encounters between nearby particles trigger chaos on a time scale 
that is typically short compared with to. As probed by A-body Lyapunov exponents computed 
in the usual way, this microchaos does not decrease with increasing A"; if anything, XN increases 
with increasing A^. However, the 'range' of the chaos, expressed relative to the total size of 
the system, does decrease with increasing A^ since the typical intcrparticle spacing scales as 
A~-^/^. Alternatively, because of the long range character of the gravitational interaction, one 
also encounters the possibility of macrochaos which will arise if the bulk density distribution 
corresponds to a bulk potential that admits global stochasticity. Both forms of chaos can play 
a nontrivial role in the A^-body problem and, as has been shown here, information about xn 
and xs can both be extracted from a judicious analysis of numerical data. 

In this sense, it would appear that although, strictly speaking, the rate at which A-body 
orbits diverge is set by the Lyapunov exponent xn, it may be misleading to assert (Heggie 
1991) that "the approximation of a smooth potential is useful for studying orbits, but not for 
studying their divergence." The smooth potential Lyapunov exponent xs does indeed provide 
useful information regarding the divergence of orbits on mesoscopic scales large compared with 
the intcrparticle spacing but small compared with the size of the system. 

Portions of this paper were written while HEK was a visitor at the Aspen Center for 
Physics, the hospitality of which is acknowledged gratefully. This work was supported in part 
by the NSF through AST0070809. 
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Figure 1: The mean spatial separation between orbits generated in frozen- A'^ systems from 
initial conditions separated in configuration space by a distance 6r{0) = 10~^. Each curve 
was generated by averaging over 100 pairs of initial conditions evolved in A^-body realisations 
of the chaotic triaxial ellipsoid plus black hole potential. The four curves correspond (from 
top to bottom) to N = 10^-^, 10^, 10^-^, and 10^. The solid line corresponds to a slope of 
0.025, generated least squares fit to the iV = 10'^ data over the interval 32 < t < 128. 
The triple dot-dashed line has a slope 0.022, equal to the mean value of the smooth potential 
Lyapunov exponent xs- The dashed line has a slope 0.75, equal to the mean value of the 
iV-body Lyapunov exponent xjv- The dot-dashed curve overlaying the data for N = 10^^^ 
represents the function Sr = A{t — to) for A = 0.008 and to = 12.0. 




Figure 2: The same as the preceding Figure, now computed for orbits in frozen- A?^ reahsations 
of the integrable Plummer density distribution. 
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Figure 3: The same data as Figure 2, now plotted on a linear scale. 
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Figure 4: The mean spatial separation between orbits with identical initial conditions evolved 
in two different frozen-AT backgrounds, each sampling the chaotic ellipsoid plus black hole 
potential. As in the preceding Figures, each curve averages over 100 pairs of orbits. The four 
curves correspond (from top to bottom) to = 10^'^, 10^, 10^'^, and 10*^. The solid line 
corresponds to a slope of 0.023, generated as a least squares fit to the N = 10^ data over the 
interval 32 < t < 160. The triple dot-dashed line has a slope 0.022, equal to the mean value of 
the smooth potential Lyapunov exponent xs- The dashed curve has a slope 0.75, equal to the 
mean value of the AT-body Lyapunov exponent xjv- 




Figure 5: The same as the preceding Figure, now computed for orbits in frozen- A?^ reahsations 
of the integrable Plummer density distribution. 




Figure 6: The mean spatial separation between orbits with identical initial conditions evolved 
in the smooth ellipsoid plus black hole potential and a frozen-AT realisation thereof. As in 
the preceding Figures, each curve averages over 100 initial conditions. The four curves again 
correspond (from top to bottom) to = lO'^'^, 10^, 10^'^, and 10^. The solid line corresponds 
to a slope of 0.023, generated as a least squares fit to the N = 10^ data over the interval 
32 < t < 160. The triple dot-dashed line has a slope 0.022, equal to the mean value of the 
smooth potential Lyapunov exponent xs- The dashed curve has a slope 0.75, equal to the 
mean value of the A/^-body Lyapunov exponent xn- 
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Figure 7: The mean spatial separation between orbits with identical initial conditions evolved 
in the smooth potential in the absence of perturbations and in the presence of energy-conserving 
white noise. The different curves correspond (from top to bottom) log^g ij = —4, —5, —6 and —7. 
The solid line has a slope 0.022, equal to the mean value of the smooth potential Lyapunov 
exponent xs- The dashed curve has a slope 0.75, equal to the mean value of the iV-body 
Lyapunov exponent xn- The dotted lines reproduce the curves in Figure 4 for AT = 10^-^ and 
10^"^, corroborating the claim that discreteness effects for N = 10?'+^/^ are well-mimicked by 
ri = 10-P. 



